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Abstract 

Diffusion of electrons in a two-dimensional system in static random magnetic 
fields is studied by solving the time- dependent Schrodinger equation numerically. The 
asymptotic behaviors of the second moment of the wave packets and the temporal 
auto-correlation function in such systems are investigated. It is shown that, in the 
region away from the band edge, the growth of the second moment of the wave 
packets turns out to be diffusive, whereas the exponents for the power-law decay 
of the temporal auto-correlation function suggest a kind of fractal structure in the 
energy spectrum and in the wave functions. The present results are consistent with 
the interpretation that the states away from the band edge region are critical. 

PACS numbers: 05.70.Fh(Phase Transitions : general aspect), 71.30 (Metal- 
Insulator transitions), 71.55.J(Localization in disordered structures), 73. 20(Electronic 
surface states) 
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§1. Introduction 

Since the development of the scaling theory of localization P],0, much work has 
been done on the transport properties of disordered electron systems |3|,[3J]. To un- 
derstand the transport properties of a random system, the concept of quantum in- 
terference plays an important role. In fact, in the last decade, many interesting and 
universal phenomena observed in mesoscopic systems have been understood in terms 
of the quantum interference [|5],||. In most cases, however, only a coherent scattering 
due to the random potentials has been discussed. Recently there has been a growing 
interest on the coherent scattering by the randomly varying magnetic fields PHlO 



In particular, considerable attention has been paid to transport properties in two- 
dimensional systems in randomly distributed magnetic fluxes. For instance, several 
authors JTT| - |r5| have discussed whether or not the Anderson transition occurs in a 
two-dimensional system in static random magnetic fields. A possibility of the An- 
derson transition in such a system has been proposed in references ||13|| , |15| and 
[|19| based on numerical calculations of conductance and the scaling analysis for the 
exponential-decay length of the Green function in quasi-one-dimensional systems. 
According to the references |13[| and [15], the mobility edge E c has been estimated as 



E c = 2.6 ~ 3.0. On the other hand, there exist numerical data which suggest the ab- 
sence of such a transition . It has been also argued quite recently that the model 



in question can be mapped to the nonlinear a model for the unitary ensemble in 
which case all states are known to be localized The conclusions of papers [|Tl~HT9| 



concerning this problem are still controversial. It is then probably fair to say that the 
question of whether the Anderson transition is realized in two-dimensional systems 
with random magnetic fluxes remains unsettled and it is a challenging problem to 
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study the transport properties of such systems. 

A two-dimensional system with inhomogeneous magnetic fields has been experi- 
mentally realized using the type II superconductors |P0|| . Theoretically, a system in 
random magnetic fields with zero mean arises in the theory of composite fermions 



that describes the fractional quantum Hall effect with half-filled Landau levels |2T 



It is therefore a very interesting problem to clarify not only the existence of the 
Anderson transition, but also the structures of the wave functions in such a system. 

In this paper we have studied dynamical aspects of the electron system in static 
random magnetic fields by solving the time-dependent Schrodinger equation numer- 
ically. This type of analysis, which is usually called the equation-of-motion method, 



has been used in the research of disordered systems P2-Z1|. In particular, the conven- 
tional Anderson model in one- and two- dimensions, having only a diagonal disorder, 



has been studied by De Raedt based on the exponential product formulae. How- 
ever, the forth-order formula used in ref. |23] involves comutators between operators 
and becuase of this it is difficult to apply it directly to the system with off-diagonal 
disorder such as a system in random magnetic fields. We therefore adopt here the 
new formula for solving the time-dependent Schrodinger equation proposed by Suzuki 



[25]. This formula has an advantage over the previous one|24]], because commutators 
between operators are not involved. It is important to note that owing to this fact 
the new formula is easily applicable to a much wider class of systems including not 
only the systems with diagonal disorder but also those with off-diagonal disorder |26| . 

The quantities we have observed are the second moment of the wave packet (r 2 (t)) c 
defined by 

(r 2 (t)) c <r 2 (t)> - (r(t)> 2 (1) 
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with 

(r n {t))= J dQr^drr^r,^ 2 , (2) 

where ip(r,t) denotes the wave function at time t, and the temporal auto-correlation 
function C(t) defined by 

C(t) = \fp{t')dt' (3) 
t Jo 

p(t) = \(m\m)\ 2 - 

Here P(t) denotes the overlap function between the initial state \ip{0)) and the state 
\ip{t)) at time t and C(t) is its time average. These quantities have been calculated 
using exact eigenvalues and eigenfunctions of the system in the cases of the Harper 
model[|27| and the critical state of the quantum Hall system [28fl. It has been shown 



that from the asymptotic behavior of the auto-correlation function C(t) we can obtain 
information about the fractal character of the energy levels and wave functions. We 
thus discuss in the present work the fractal character of the wave function, which 
is usually observed just at the critical point or at length-scales much smaller than 
the characteristic length of the system. Analyses of the fractal character of the wave 
functions at the critical energy of the Anderson transition have been made for two- 
dimensional systems in a strong magnetic field |f29| - [32]j and also for three-dimensional 



systems |33|-p6| . In particular, the relationship between transport properties and the 



fractal structures of the system is a very interesting problem |37| . 

The organization of the present paper is as follows. In the following section, we 
explain our method which is based on the higher-order decomposition of exponential 
operators. The expected asymptotic behaviors of the second moment of the wave 
packet and the auto-correlation function are also discussed in this section. In section 
3, the above method is applied to the two-dimensional systems in static random 



magnetic fields and the numerical results for the growth of the second moment and the 
decay of the auto-correlation are presented. The final section is devoted to discussions 
of the physical implication of our numerical results. 



§2. Numerical Method 



In order to solve numerically the time-dependent Schrodinger equations, we have 



adopted the method based on the higher-order decomposition of exponential operators p5 
As is seen below, the method is quite general and easily applicable to complicated 
systems described by the tight-binding Hamiltonian with inhomogeneous hopping 
elements. The basic formula we have used is the fourth-order decomposition of ex- 
ponential operators: 

exp[x(A 1 + A 2 + • • • + Ag)] = S(xp)S(x(l - 2p))S{xp) + 0(x 5 ), (4) 

where 

and the parameter p is given by p = (2 — v^) -1 . Here Ai,...,A q are arbitrary 
operators. 

Let us consider, in general, the tight-binding Hamiltonian on the two-dimensional 
square lattice that has nearest neighbor hopping amplitudes {tij;t*j = tj t i} and the 
diagonal elements {ei}: 

H = ~ E • Y, :(';(:■ (6) 

<i,j> i 

where C\{Ci) denotes a creation (annihilation) operator of an electron at the site i. 
We then decompose this Hamiltonian into five parts (see figure 1), namely, 

H = £#n, (7) 



r^Godd r s 

= ^ tf+xjrCjt+.tfi'r + h.C. , 

r x £even r y 

H3 = tf + y^ci + ^Cf + h.c 

^ r a £odd 
T"sc r H Seven 

# 5 = Y^etClC?. (8) 
f 

In equations (^|), jc(y) denotes the unit vector in the x(y)— direction and r x (r y ) G Z 
the 37(2/ ) component of the position vector r of the sites. All the length-scales are 
measured in units of the lattice constant do- The Hamiltonian Hi, for instance, 
describes the hoppings between the sites whose positions are r and r + x for r x = 
1, 3, 5, . . . and r y = 1, 2, 3, . . .. Note that each Hamiltonian H n consists of operators 
which commute with each other, and hence we can obtain the analytical expressions 
for exp(—iH n t/H) by diagonalizing two by two matrices. According to the formula 
(H|), it is easy to see that the state vector \tp(t + 8t)) at time t + St can be obtained as 

\^(t + 5t)) = exv(-i8tH/h)\ip(t)) 

= S 2 (-i6tp/h)S 2 (-i6t(l-2p)/h)S 2 (-i5tp/h)\t(;(t)) 

+0{5t 5 ) (9) 

with 

S 2 ( x ) = e ^/2 e xH 2 /2 . . . e xHt/2 e xH Be zHt/2 . . . ^2^/2 (1Q) 

It should be noted that in the present approach the unitarity of the process is obvi- 
ously preserved and that we have high enough accuracy to perform large time-scale 
simulations. It is important to note that in this method we do not need to diago- 
nalize the whole Hamiltonian. We can thus simulate larger systems, for example a 
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499 by 499 square lattice, with high precision, which is crucial if we are to extract 
meaningful properties from the numerical data for random systems. 

In order to consider the time evolution of the wave packet with fixed energy E, 
we have carried out a numerical diagonalization of a subsystem H^ whose size is N 
by iVo, located at the center of the whole system and have chosen an eigenstate of 
i/jVp with eigenvalue En ~ E as the initial wave function. 

Next, let us consider the asymptotic behaviors of the auto-correlation function 
C(t) and the second moment of the wave packet (r 2 ) c in several regions, namely, the 
metallic, the localized (insulating) and the critical regions. First, in the localized 
region it is clear that the second moment remains finite (r 2 (t)) c — » d(d + l)£ 2 /4 
in the limit as t — > oo, where £ denotes the localization length that describes 
the exponential-decay of the wave function as IV'COI ^ ex P[ — r /£]- Obviously, the 
quantity C(t) also goes to a nonzero value in the same limit. In contrast, in the 
metallic region it is expected that the second moment grows in proportion to time t: 

(r 2 ) c = 2dDt (11) 

and that the auto-correlation C(t) decreases as 

C(t) oc r d/2 . (12) 

Here, the coefficient D denotes the diffusion coefficient and d the dimensionality of 
the system. 

In the critical region where the wave function has a fractal structure, a non- 



conventional behavior of the auto-correlation function C{t) has been observed fl27|,|28 



In such a case, the auto-correlation function decays in proportion to t a for t > 1 



with a ^ d/2 reflecting the fractal character of the energy spectrum |27j and of 
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the wave function fl28|. It has been shown that the exponent a is related to the 



generalized dimension D 2 of the energy spectrum by a = D 2 [57], and, further, that 
in the case of the quantum Hall system it is also related to the generalized fractal 



dimension D 2 [|3T| of the wave function by a = D 2 /2 pB |. The generalized dimension 
D 2 of the energy spectrum is defined as the scaling exponent for the correlation of 
energy levels: 

d/j,^(Q)(u) d/z^( )(u/) oc l D2 as I -> 0, (13) 

J w—l 

where iim$\{oj) denotes the spectral measure with respect to an initial state ip(0) 



27J. The generalized fractal dimension D 2 of the wave function is, on the other 



hand, defined by: 
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P 2 (A,£) = £ E mr,E)\ 2 ) ocA D2 as A - 0, (14) 
i y ren t (\) J 

where {f2j(A)} denote AA^ x AA^ boxes that cover the whole system (N x N) p8| . 
On the basis of these relations, the non-conventional behavior of the auto-correlation 
function C (t) can be considered as evidence for fractal structure in the system. Since 
the multifractal behavior has been demonstrated not only at the critical point, but 
also for length-scales smaller than the characteristic length, such as the localization 
length in the localized region fl35| , a similar type of non-conventional behavior of C(t) 



is expected if we consider length-scales much smaller than the characteristic length 
of the system. 
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§3. Model and Numerical Results 

The two-dimensional system with random magnetic fluxes is described by the 
tight-binding Hamiltonian: 

H = ~ E VijCjCj, (15) 

<i,3> 

where 

V iti = Ve Xp [i0 iJ ]. (16) 
The phases {9ij} are related to the magnetic fluxes {(pi} through the relation: 

@i+x,i + @i+x+y,i+x + @i+y,i+x+y + &i,i+y = ~ 2ll§i/ fa, (17) 

where fa and fa = hc/\e\ denote the magnetic flux through the plaquette + 
x,i + x + y,i + y) (figure 2) and the unit flux, respectively. It should be noted here 
that we have the gauge degrees of freedom in determining the phases {9ij} for given 
fluxes {fa}. The magnetic fluxes through each plaquettes of the square lattice are 
distributed uniformly in the interval [— fa/2, fa/2], and the fluxes are assumed to 
have no spatial correlation, i.e. 

(0i0j)av = j^ttj- ( 18 ) 

Here the angular bracket ( ) av denotes the average over the distribution of fluxes. 

We have calculated the auto-correlation function C(t) and the second moment of 
the wave packet (r 2 ) c at various energies. The sizes of the systems are 499 by 499 
for the energies near the band center and 199 by 199 for the energy near the band 
edge. We have carried out an exact diagonalization for the 21 by 21 subsystem at the 
center of the system and taken the eigenfunction of the subsystem whose eigenvalue 
is closest to the given energy E as the initial wave packet. By this procedure we can 
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simulate the diffusion of the wave packet whose energy is approximately equal to E. 
The single time step St is taken to be 0.2(h/V) in the simulation. With this choice 
of St, fluctuations of the expectation value of the Hamiltonian 

(h) = (mmm) (19) 



can be safely neglected p4j throughout our simulations (t < 8000(7i/V)). In the 
actual simulation, the quantities (r 2 ) c and log C(t) are averaged over five realizations 
of random magnetic flux distribution. 

Let us first discuss the diffusion near the band edge. The states near the band edge 



are expected to be localized Hllip^| . Our results for the energy (H) ~ — 3.35V are 



shown in figures 3(a) and 3(b), which is based on a calculation of five configurations of 
random magnetic fields. The second moment of the wave packet {r 2 ) c seems to have 
a finite value in the limit t — > oo. The localization length can be roughly estimated 
as £ ~ 14, which is consistent with the value obtained from the finite-size scaling 
analysis of ref. [14]. It may not be clear from the behavior of the second moment of 



the wave packet alone whether the states are indeed localized or not. The behavior 
of the auto-correlation C(t), however, clearly suggests that the state at time t has a 
nonzero overlap integral with the initial state in the limit t — > oo and thus the states 
are localized in this energy region. 

Next, we discuss the diffusion near the band center. The results for the energies 
(H) /V ~ 0, — 1, —2, —2.6 are presented in figures 4, 5, 6 and 7, respectively. The 
growth of the second moment of the wave packets turns out to be diffusive; the 
numerical value of the exponent a for the power-law growth defined in 

(r 2 (t)) c oc t a (20) 

is very close to 1. It should be noted here that, in the critical case of the Harper 
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model, the exponent a has been estimated to be about 0.97 [38|, 59"|, which is also close 
to 1 [4"0|. The exponent a has been also estimated to be 1 for the critical states of the 
quantum Hall system p8 |. By assuming a linear dependence of (r 2 (t)) with respect 



to t, we can evaluate the diffusion coefficient D defined in eq.(|TT]), which is related to 
the conductance of the system. The obtained values of the diffusion coefficients for 
energies —2.6, —2, —1 and are consistent with those obtained using the Landauer 



formula 44 



On the other hand, the exponents a for the power-law decay of the auto-correlation 
function show deviations from the value d/2 = 1 in the same energy region. For the 
energy region —2.6 < E < 0, the exponent a is estimated to vary from 0.9 to 0.88, 
which is smaller than 1. These deviations suggest the presence of fractal properties 
of the energy spectrum and the wave function, just as observed in the cases of the 
Harper model and the quantum Hall system P7|, Z8|. The deviations themselves are 



smaller than that obtained for the critical case of the quantum Hall system p8 
suggesting a larger fractal dimensionality of the present system. 

§4. Discussions 

We have performed numerical simulations for the diffusion of electrons in static 
random magnetic fields using a new formula for the decomposition of exponential 



operators 25 . The systems we have treated are much larger than those in previous 



studies of this kind [p4], pq|. The present method reproduces the results previously 
found from the finite-size scaling analyses, namely that the states near the band edge 
are localized. On the other hand, away from the band edge, the second moment 
of the wave packets increases linearly in time, which suggests that the conductance 
is non-zero in this region. If the system is described by the conventional diffusion 
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equation, the exponent a should be d/2 = 1 and in that case no fractal structure 
is present. However, the obtained values of the exponents a of the power- law decay 
of the auto-correlation function are smaller than d/2 — 1. It is expected from these 
deviations of the values of a from 1 that some kind of fractal structure exists in the 
energy spectrum as well as in the eigenfunctions of the present system. This type 
of behavior, namely the linear increase of (r 2 ) c in time and the non- conventional 
exponents for the power-law decay of the auto-correlation function, has been observed 
in the critical case of the Harper model and the critical states in the quantum Hall 



system [27],|28]. It may therefore be natural to suppose that the states away from the 
band edge in the present system with random magnetic fluxes are also critical. 

It has been shown that, by assuming conformal invariance, the exponent rj for 
the power-law decay of the critical wave function in two-dimensional systems can be 
related to the exponential decay-length of the Green function in quasi-one dimensional 
strip systems PHfI2|. The exponent r], which describes the power-law decay of the 
density correlation by (|-?/'(r)| 2 |'?/'(0)| 2 ) av ~ \r\~ 2r > in two dimensions, is expressed as 

m 



* = ik- (21) 

where A c denotes the renormalized localization length, defined by 

A c = Km £ M (E)/M. (22) 

M—*oo 

Here £m{E) denotes the localization length along a strip of width M with energy 
E corresponding to the critical states in infinite two-dimensional systems [43[|. The 



fractal dimensionality d* of the wave function is then obtained from the exponent rj 

m 

d* = d- 2r]. (23) 



as 
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In this way, the fractal dimensionality d* of the wave function in a two-dimensional 
system at its critical energy, or at the length-scales much smaller than the characteris- 
tic length, can thus be evaluated from the scaling behavior of the localization length 
along the quasi-one-dimensional systems. Making use of this theory we have 
performed independently the calculation of £m for a few random field configurations 
in strips whose sizes are M x L with L = 10 5 and M = 8, 12, 16, 20, 34, 40 and 50 in 
units of the lattice spacing a (figure 8) and have estimated the exponent r] through 
the formulae (pi]) and fl22|) . The values of d* obtained by the formula (|23|) are then 
compared with the generalized fractal dimension D 2 obtained from the exponent a 
by assuming the relation D 2 = 2a PS| (Table 1). It is rather remarkable that the 
values of d* and D 2 coincide fairly well with each other. 

We have thus argued, based on our numerical results, that there is evidence for 
fractal structure in the energy spectrum and the eigenfunction in the region away 
from the band edge. However, we admit that the fractal character of the present 
system turns out to be not so strong compared with the case of the critical states in 



the quantum Hall system |[28|| . At the band center, in particular, we cannot exclude 
the possibility that there is no fractal structure, because the value of the exponent 
a turns out to be rather close to 1. Since the fractal structure is observed when 
the characteristic length is much larger than the size of the system, then precisely 
speaking, our results do not distinguish the critical states from the localized states 
with extremely large localization length. Simulations in much larger systems and 
the investigation of the probability distributions of (r 2 ) c and C(t) will be needed in 
order to obtain more detailed information about the thermodynamic properties of 
this system. 

Finally, we emphasize that the formula we have used in this work is quite general. 
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It would be very interesting to study many other dynamical problems, such as the 
anomalous diffusion [OJ and the mesoscopic dynamical echo [[16j by making use of 
the present method. 
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E 


A c 


d* = d - 2r] 


D 2 = 2a 





3.04 ±0.04 


1.791 ±0.002 


1.8 ±0.06 


-1 


3.00 ±0.08 


1.788 ±0.006 


1.8 ±0.06 


-2 


2.73 ±0.07 


1.766 ±0.006 


1.8 ±0.03 


-2.6 


2.34 ±0.09 


1.73 ±0.01 


1.76 ±0.02 



Table 1: Fractal dimensions obtained from the exponents r] and those from the ex- 
ponents a for E = 0, —1, —2, —2.6. The values A c are estimated using M x L strips 
with 8 < M < 50 and L = 10 5 . 
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Figure Captions 



Fig. 1: The decomposition of the Hamiltonian. The thin solid, thick solid, thin 
dashed and thick dashed lines correspond to the Hamiltonian H\, H2, H3 
and if 4 in equation (|J), respectively. 

Fig. 2: The plaquette (i, % + x, i + x + y, % + y) and the flux & through it. 

Fig. 3: (a) The growth of the second moment (r 2 (t)) c of the wave packet and (b) 
the time-dependence of the auto-correlation function C (t) for the energy 
(H) ~ — 3.35(— 3.35 ± 0.01)K. The bars around the data points in figures 
3, 4, 5, 6 and 7 indicate the fluctuations of the mean values with respect 
to five realizations of random magnetic field configurations. 

Fig. 4: (a) The growth of the second moment {r 2 (t)) c of the wave packet and 
(b) the time-dependence of the auto-correlation function C(t) for the 
energy (H) ~ 0(±1.0 x 10~ A )V. The solid line in (a) corresponds to 
D = 1.59alV/h. The exponent a is estimated to be 0.90 ± 0.03. For 
comparison, the conventional decay of C(t) oc t~ x is also shown by the 
dashed lines in figures 4(b), 5(b), 6(b) and 7(b). 

Fig. 5: (a) The growth of the second moment (r 2 (t)) c of the wave packet and (b) 
the time-dependence of the auto-correlation function C(t) for the energy 
(H) ~ — 1(— 1.01 ± 0.02)V. The solid line in (a) corresponds to D = 
1.54agV /h. The exponent a is estimated to be 0.90 ± 0.03. 

Fig. 6: (a) The growth of the second moment {r 2 (t)) c of the wave packet and (b) 
the time-dependence of the auto-correlation function C (t) for the energy 
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(H) ~ -2(-1.99 ± 0.01)V. The solid line in (a) corresponds to D = 
1.33alV/h. The exponent a is estimated to be 0.90 ± 0.014. 

Fig. 7: (a) The growth of the second moment {r 2 (t)) c of the wave packet and (b) 
the time-dependence of the auto-correlation function C (t) for the energy 
(H) ~ — 2.6(— 2.61 ± 0.01)V. The solid line in (a) corresponds to D — 
0.98alV/h. The exponent a is estimated to be 0.88 ± 0.01. 

Fig. 8: Energy dependence of the exponential decay length £m divided by the 
width of the strip M |L4|,[|3||. The circle, square, diamond, triangle, bullet, 
painted square and painted diamond correspond to M — 8.12, 16, 20, 30, 40 
and 50 respectively. The length of the strip L is 10 5 and the average is 
performed over more than two random flux configurations. Due to the 
exact particle-hole symmetry, only the negative energy region is plotted. 
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